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ABSTRACT 


Recursion  equations  are  developed  for  solution  of  transient  tem¬ 
perature  distributions  in  an  infinite  hollow  cylinder  with  nonsteady 
boundary  conditions.  The  solution  is  shown  to  be  applicable  to  any 
imposed  boundary  condition  and  is  also  shown  to  be  able  to  handle  the 
special  case  of  the  solid  cylinder.  A  computer  program  is  written 
and  applied  to  two  examples.  A  comparison  of  the  numerical  results 
with  classical  exact  solutions  reveals  close  agreement  between  the 
two  types  of  solutions  for  the  particular  cases  considered. 
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SECTION  I 
INTRODUCTION 


The  design  of  industrial  processes  frequently  requires  a  knowledge 
of  the  transient  heat -transfer  phenomena  within  infinite  cylinders.  The 
heat  treatment  of  solid,  long  steel  billets,  or  the  analysis  of  thermal 
stresses  in  pipelines  carrying  hot  or  cold  fluids  and  subject  to  nonsteady 
flow  conditions  are  examples  where  an  analysis  of  transient  heat  trans¬ 
fer  is  desired. 

The  derivation  of  the  governing  differential  equation  for  transient 
heat  flow  in  a  cylinder  is  well  known  and  cam  be  found  in  texts  such  as 
Schneider  (Ref.  1)  or  Eckert  and  Drake  (Ref.  2).  In  the  most  general 
case,  the  cylinder  is  hollow,  has  an  arbitrary  internal  distribution  of 
temperature,  and  is  subject  to  nonsteady  environments,  not  necessarily 
the  same,  at  both  the  inside  and  outside  surfaces.  Appendix  I  contains 
the  derivation  of  the  governing  differential  equation,  together  with  a 
formal  statement  of  this  general  set  of  boundary  conditions. 

Solutions  to  the  general  equation  subject  to  a  more  restrictive  set 
of  boundary  conditions  have  been  developed  by  several  authors.  Schneider 
(Ref.  1)  develops  solutions  for  solid  cylinders  with  arbitrarily  specified 
initial  radial  temperature  distribution  and  with  a  temperature,  T0,  sud¬ 
denly  applied  to  the  outside  face.  Schneider  also  develops  solutions  for 
the  solid  cylinder  subject  to  a  uniform  initial  temperature  and  suddenly 
immersed  in  a  fluid  at  temperature  Tg.  Boelter,  Cherry,  et  al.  (Ref.  3) 
also  present  solutions  for  transient  temperature  distributions  in  solid 
cylinders,  with  arbitrary  initial  temperature  distributions,  that  are  sud¬ 
denly  immersed  in  a  constant  temperature  fluid.  However,  analytical 
solutions  of  the  differential  equation  subject  to  the  boundary  conditions  of 
Appendix  I  do  not  appear  to  be  available  in  the  literature. 

This  document  addresses  itself  to  the  problem  of  developing  a 
numerical  technique  for  solving  the  general  differential  equation  subject 
to  the  general  set  of  boundary  conditions. 


SECTION  II 

DERIVATION  OF  THE  RECURSION  EQUATIONS 


The  recursive  technique  is  introduced  by  dividing  the  cylinder  into 
N  cylindrical  shells,  each  of  equal  thickness,  2AR.  Next  the  assumption 
is  made  that,  for  sufficiently  small  increments  of  time,  the  problem 
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approaches  a  quasi-steady  state  process  for  each  of  the  incremental 
elements.  Consequently,  by  making  successive  energy  balances  on 
small  elements  of  the  pipe  for  small  time  periods,  a  steady-state  ap¬ 
proximation  of  the  actual  process  can  be  made. 

Figure  1  (Appendix  II)  shows  a  unit  length  of  cylinder  divided  into 
N  concentric  elements  of  equal  thickness.  Choosing  N  large  enough 
allows  each  elemental  shell  to  be  considered  isothermal.  If  the  cylinder 
is  considered  a  composite  of  N  isothermal  shells,  each  at  its  own 
specific  temperature,  then  over  short  periods  of  time  the  temperature 
change  of  any  element  will  be  a  function  only  of  the  temperature  of  that 
element  and  the  temperature  of  adjacent  elements. 

Referring  to  Fig.  1,  the  geometrical  aspects  of  the  pipe  indicate 
that  the  thickness,  t,  of  each  shell  is  given  by  Eq.  (1) 

t  =  (R0  -  Ri)  +  N  =-  2  A  R  (1) 

where  Ri  and  R0  are  the  inside  and  outside  pipe  radii,  respectively,  and 
N  is  the  chosen  number  of  concentric  shells. 

An  energy  balance  is  now  made  on  the  first  cylindrical  shell: 

Rate  of  Energy  In  (from  Gas)  -  Rate  of  Energy  Out  (to  Element  2)  =  Time  Rate  of 
Change  in  Energy  Stored  in  Element  1. 


Expressed  mathematically  the  energy  balance  becomes 
l.  /a)  T  1  ~  Tj)  r  u  dT,  r  r  V  ^ 1 

hj(0) A1(Tglw)  -Tl)  +  Hj/Rj  -  ~  CpV‘  — IT~ 


hi (0)  2nRj  (Tgi(0)  -  T.) 


2  n  K(Tj  -  T.) 
fn  Rj/R, 


=  C  p  n  [(Ri  +  2AR)J  -  R;2] 


(T'  -  T.) 

W~ 


(2) 


Here  Tj  is  taken  to  represent  the  new  temperature  of  Element  1 

after  exchanging  energy  with  adjacent  elements  during  the  time  period 
A0.  C  is  the  pipe  metal  specific  heat,  K  is  the  pipe  metal  thermal  con¬ 
ductivity,  p  is  the  pipe  metal  density,  Vj  is  the  volume  per  unit  length 
of  the  first  (inside)  cylindrical  shell,  and  Ai  is  the  inside  area  of  the 
cylinder  per  unit  length. 


Recalling  the  definition  of  AR  allows  Rj  and  R2  to  be  expressed  as 
follows : 


Rt  =  Ri  +  A  R  =  Ri  +  (R0  -  Ri)  /2N 


R2  =  Rj  +  3  A  R  =  Rj  + 


3(R0  -  R;) 
2N 


(3) 
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Substituting  Eqs.  (1)  and  (3)  into  Eq.  (2)  gives: 


jhi(0)  2n-Ri  (Tgi(0)  -  Tj)  +  2  n  K  (T2  -  T,)/£n 


=  C 


-  r^ct;  -  T,) 


\ 


Expanding  the  squared  terms  and  simplifying  yields 
I2  rhj(0)Ri  (Tgi(0) 


t.)  ♦  2.Kcr,  -  T.yf » [;" ;*•][ « 


=  c 


p  n 


[R02  -  (2N  -  2)  R„  R,  +  (1  -  2N)  Rj2]  (Tj  -  T,) 
__ 


(4) 


It  now  becomes  convenient  to  redefine  certain  parameters  into 
groups.  Let  J3  =  R0/K,  let  r  represent  the  ratio  of  inside  radius  to  out¬ 
side  radius  (Ri/R0),  and  let  4>  be  as  defined  below: 

,  2  K  N2  A0  2  N*  A0  a 

ijj  =  -  »  - 

C  p  R02  R„2 

where  a  is  the  thermal  diffusivity. 


Substituting  these  parameters  into  Eq.  (4)  and  solving  for  the  new 
element  temperature,  T^,  gives 


TJ  =  T,  +  1I1 


jr|3  hj(0)  [Tgj(0)  —  Tj]  4.  (T2  -  T,)/fn 

'(2N  -  3)  r  +  3]  'i 

(2N  -  1)  r  +  1  ( 

)  (]  -  r)  [(2N  -  1)  r  +  i 

1 

(5) 


Looking  now  at  the  inner  elements  of  the  cylinder  and  writing  an 
energy  balance  for  Element  2  gives 


2 v  K  (T,  -  Tj)  277  K  <Tj  -  T2) 
+ 


tn  (R2/R,) 


in  (Rj/R2) 


a  e  =  c  P  v2  (t;  -  t2) 


(6) 


Generalizing  Eq.  (6)  to  the  kth  cylindrical  element  gives 

A0  =  CpVk  (T'k-Tt) 


277  K  (Tk_!  -  Tk)  2  77  K  {Tk+1  -  Tk) 
- - -  + 


^(Rt/Rk.j)  fn(Rfc+1/Rk) 


(7) 
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Solving  for  gives 


T'k  =  Tk  + 


2rrA0K 
Cp  Vk 


Tk-  1  ~  Tk 

£n  (Rk/Rk_  !> 


Tk  -  i  ~  Tk 
tn  (Rk  +  i/Rk) 


(8) 


Substitution  for  V^,  In  (R^/R^-i),  lnfR^i/R^),  and  applying  the 
definition  of  V>  gives 


Tk  =  Tk  +  <jf 


/  Tk_!-Tk 

Tk  +  I-  Tk  \ 

\  ^  ["  (2k-  1)  -  (2N  -  2k  +  l)r" 

tn 

"(2k  1>  +  (2N  -  2k  -  l)r“|  I 

)  |_  <2k  -  3)  +  (2N  -  2k  +  3>r 

_(2k  -  1)  +  (2N  -  2k  +  l)rj  \ 

\  (1 -r)  [(2N  -  2k  +  l)r+ (2k- D]  ^ 

(9) 


Next  a  look  is  taken  at  the  outermost  cylindrical  element.  The 
energy  balance  gives 

/  \ 

-INo  ~  JW)  +  2tt  h0(fl)  R0  (T go(0)  -  TN)(  A  0  =  C  p  VN(Tn  '  -  TN)  ( 10) 

I  ^  _Ra-3AR_  ) 


|E  TRo  (2N  -  l)l+  Ri  1  +  2,rh"(fi)R»  1Ti»(0)  ~  TN}|  A6  = 
LRo  (2N  -  3)  +  3  R;  J  ) 


C  p  rr |ros  -  |~R0  -  HV-  TN)|  ^ 


1) 


Solving  for  gives 


Tn  =  Tn  +  ^ 


tn_ i-tx 


2N  -  l  +  t 


2  N  **-  3  +  3  r  I 


0ho{0)[Tgo(0)  -  TN] 


(1  -  r)  [r+  2N  -  l] 


(12) 


Equations  (5),  (9),  and  (12)  will  yield  the  new  temperature  of  the 
inner  element,  the  middle  elements,  and  the  outer  element,  respec¬ 
tively,  after  a  given  time  interval  A0.  The  information  required  to 
solve  for  the  new  element  temperatures  is  as  follows: 
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1.  The  previous  element  temperature  as  well  as  previous  tem¬ 
peratures  of  adjacent  elements. 

2.  Time -dependent  functions  for  inside  and  outside  surface  film 
coefficients  (if  applicable). 

3.  A  specification  of  the  outside  radius  of  the  cylinder  and  the 
radius  ratio. 

4.  A  specification  of  the  material  properties  of  the  cylinder. 


SECTION  III 

CONVERGENCE  AND  STABILITY 


.3.1  CONVERGENCE 


In  a  qualitative  sense,  convergence  is  a  measure  of  the  ability  of 
the  numerical  approach  to  yield  results  which  agree  with  those  results 
that  could  hypothetically  be  determined  from  an  exact  solution.  If  in 
Eqs.  (5),  (9),  and  (12)  the  value  of  A0  was  allowed  to  approach  zero  and 
the  value  of  N  was  allowed  to  approach  infinity,  then  the  numerical  solu¬ 
tion  would  converge  to  exactly  the  same  results  as  the  exact  solution. 
However,  it  is  impossible  to  allow  the  limiting  cases  for  A0  ,and  N  to 
exist  or  even  to  be  realistically  approached.  Hence,  convergence  is  the 
measure  of  how  close  the  exact  solution  can  be  approached  using  finite 
values  of  A0  and  N.  In  the  interest  of  time,  it  behooves  us  to  use  the 
largest  value  of  A0  and  the  smallest  value  of  N  that  will  give  convergent 
results . 


Convergence  is  difficult  to  study  analytically.  However,  a  simpli¬ 
fied  approach  can  lead  to  criteria  for  establishing  limits  to  the  con¬ 
stants  A 6  and  N  beyond  which  the  solution  will  begin  to  diverge. 

Equation  (5)  may  be  rewritten  as  follows: 


t;  =  t, 


-  * 


hi(0)  + 


1 


"( 2  N  —  3 !  r  -  3~ 


<2M  -  l)r+  1 


(1  -  r)  L(2N  -  1)  r+  1|] 


l  r/3hi(0)Tgi(0)  + 


Ir 


+  */ _ _ _ UjN  - 

V  \  (1  -  r)  l(2N  -  1)  r+  l] 


(2N  -  3)r  +  3| 


(13) 
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Inspection  of  Eq.  (13)  indicates  that  it  can  be  expressed  in  the  following 
form: 

T'x  =  Tt(l  -  At)  +  B  (14) 


where  the  terms  Ai  and  B  are  complicated  terms  involving  the  param¬ 
eters  r,  j3,  h^e),  Tg£(0),  T2,  and  N.  It  is  evident  that  the  term  B 

in  Eq.  (14)  will  always  be  positive  while  the  term  (1  -  Ai)  can  be  either 
positive  or  negative  depending  on  the  value  of  Aj_.  If  A  ^  is  greater  than 
unity,  then  (1  -  Aj)  becomes  negative,  and  the  solution  diverges.  Solu¬ 
tion  divergence  occurs  because  for  (1  -  Ai)  <  0,  with  the  value  of  B 
fixed,  any  new  estimate  of  Ti  is  inversely  related  to  T^.  Hence,  higher 
values  of  Ti  call  for  lower  values  of  Ti  and  vice  versa,  which  is  a 
violation  of  intuitive  reasoning  and  in  some  instances  could  entail  a 
violation  of  the  second  law  of  thermodynamics.  Thus,  for  the  first  in¬ 
cremental  element,  a  convergence  criterion  would  restrict  values  of  A^ 
toAi|  1.  Equation  (15)  is  the  formal  statement  of  convergence  criterion 
based  upon  the  first  element. 

rj8hi<0)  +  — 
in 

(1  -  r)  [(2N  -  l)r+  l]  (’ 


[(2N  -  3)r  +  3~|  i 
(2N  -  l)r  +  lj  f 


<  1 


Following  similar  reasoning,  convergence  criteria  can  be  devel¬ 
oped  for  the  kth  and  Nth  elements,  Equations  (16)  and  ( 17)  are  formal 
statements  of  the  convergence  criteria  for  the  kth  and  Nth  elements, 
respectively. 


Ak  .=  t 


I 

i 

"(2k  -  1)  +  (2N  -  2k  +  l)r" 

LlH 

"(2k  +  1)  +  (2N  -  2k  -  l)r 

p  1 

_(2k  -  3)  +  (2N  -  2k  +  3)r_ 

_(2k  -  1)  +  (2N  -  2k  +  l)r_ 

(1  -  r)  [<2N  -  2k  +  l)r  +  (2k  -  1}] 


<  1 


(16) 


An 


"2-V  +  r  -  1" 

_2N  -  3  r  -  3. 

(J  -  r)  (2N  4 

i 

/ 

<  1 


(17) 


Figure  2  shows  a  typical  application  of  the  convergence  criteria. 

The  recursion  equations  of  Section  II  were  applied  to  a  situation  where 
a  solid  cylinder  was  suddenly  immersed  in  a  fluid  at  a  constant  tempera¬ 
ture  Tg.  The  film  coefficient  was  held  constant.  For  this  particular 
application,  the  inside  radius  in  the  recursion  equations,  along  with  the 
inside  film  coefficient,  were  simply  set  equal  to  zero.  Plotted  on  Fig.  2 
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is  the  exact  solution  of  the  problem  along  with  the  results  of  two  recur¬ 
sive  solutions  —  one  in  which  the  convergence  criteria  is  satisfied  and 
one  in  which  the  convergence  criteria  is  not  satisfied.  Inspection  of 
Fig.  2  indicates  that  the  recursive  solution  in  which  the  convergence 
criteria  is  satisfied  converges  to  the  exact  solution.  However,  for  the 
case  of  nonsatisfaction  of  the  convergence  criteria,  the  recursive  solu¬ 
tion  is  highly  divergent. 


3.2  STABILITY 

Stability  is  a  measure  of  the  ability  of  the  numerical  system  to  ab¬ 
sorb  systematic  errors  caused  by  roundoff  and  by  the  finite  element 
assumption.  A  stable  system  damps  out  errors  and  distributes  them 
evenly  throughout  the  system.  Like  convergence,  stability  is  not  easy 
to  analyze.  However,  Smith  {Ref.  4)  gives  a  well-developed  analysis 
of  stability  for  certain  classes  of  numerical  problems.  Smith  indicates 
that  one  method  of  gaining  insight  into  the  stability  of  a  numerical  solu¬ 
tion  is  to  apply  the  numerical  equations  to  the  errors  themselves. 

Hence,  by  assuming,  for  example,  a  unit  error  somewhere  in  the  sys¬ 
tem  and  watching  the  propagation  of  this  error  as  the  solution  progresses, 
an  indication  of  system  stability  becomes  available.  Table  I  (Appen¬ 
dix  IV)  is  a  tabulation  of  the  propagation  of  a  unit  error  that  is  assumed 
to  occur  at  time  equal  zero,  in  the  tenth  element,  of  a  twenty-element 
system. 

Table  I  reveals  that  errors  in  the  system  do  indeed  tend  to  distribute 
themselves  and  to  diminish  with  time.  Although  a  check  of  this  sort  is 
not  an  all-conclusive  proof  of  stability,  it  does  indicate  that  the  system 
tends  to  be  stable,  provided  roundoff  error  is  not  excessive  and  can  be 
expected  to  be  random  with  regard  to  sign. 


SECTION  IV 
COMPUTER  PROGRAM 


Equations  (5),  (9),  and  (12)  were  programmed  for  solution  on  the 
IBM  360/50  computer.  Programming  was  carried  out  in  G  Level 
FORTRAN. 

The  program  is  set  up  to  take  any  type  of  boundary  condition  or 
initial  condition.  The  program  is  expected  to  be  called  on  to  solve  the 
problem  for  two  general  types  of  boundary  conditions: 
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1.  Type  I  boundary  conditions  involve  an  arbitrary  specification 
of  initial  temperature  distribution  within  the  cylinder  wall, 
together  with  the  specification  of  a  suddenly  applied  tempera¬ 
ture  to  either  or  both  the  inside  and  outside  cylindrical  sur¬ 
faces.  The  initial  temperature  distribution  is  input  as  data, 
while  the  suddenly  applied  surface  temperature  is  included  in 
the  program  as  an  arithmetic  statement  function.  The  sud¬ 
denly  applied  surface  temperature  can,  in  the  general  case, 
be  time  variant. 

2.  Type  II  boundary  conditions  also  involve  an  arbitrary  specifi¬ 
cation  of  initial  temperature  distribution  within  the  cylinder. 
However,  Type  II  problems  differ  from  Type  I  problems  in 
that  the  boundary  conditions  are  specified  in  terms  of  a  time- 
varying  surface  film  coefficient  and  fluid  temperatures  adja¬ 
cent  to  each  surface  of  the  cylinder.  Both  the  inside  and  out¬ 
side  surface  film  coefficients  are  included  in  the  program  as 
time-dependent  arithmetic  statement  functions.  Likewise,  the 
inside  and  outside  fluid  temperatures  are  included  as  time- 
dependent  arithmetic  statement  functions.  Film  coefficients 
and  fluid  temperatures  can,  as  a  special  case,  be  considered 
constant.  The  surface  film  coefficients  can  also  be  considered 
as  functions  of  other  parameters  such  as  Quid  temperature 
and  the  thermodynamic  and  transport  properties  of  the  fluid,  if 
desirable.  The  thermodynamic  and  transport  properties  of  the 
fluid  must  then  be  specified  and  included  in  the  program  as 
functions  of  time  or  temperature. 

A  combination  of  Type  I  and  Type  II  boundary  conditions  can  also 
be  handled  by  the  program. 

The  special  case  of  the  solid  cylinder  is  handled  readily  by  the 
program  if  the  radius  ratio  is  simply  set  equal  to  zero. 

Appendix  III  gives  a  FORTRAN  listing  of  the  program  as  applied  to 
Sample  Problem  Two  of  the  following  section. 


SECTION  V 

APPLICATION  OF  THE  RECURSIVE  METHOD 


The  application  of  the  numerical  method  will  be  illustrated  by  two 
sample  problems.  Sample  Problem  One  involves  transient  heat  transfer 
in  a  solid  cylinder  with  Type  II  boundary  conditions  at  the  outside  surface. 
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The  results  of  the  recursive  solution  are  compared  with  exact  solutions 
from  the  literature  for  the  same  problem. 

Sample  Problem  Two  is  a  typical  example  of  the  application  of  the 
method  to  a  hollow  cylinder  with  time -dependent.  Type  II,  boundary  con¬ 
ditions  on  both  the  inside  and  outside  surfaces. 


5.1  SAMPLE  PROBLEM  ONE 

Schneider  (Ref.  1)  develops  the  exact  solution  for  the  transient  tem¬ 
perature  distribution  in  a  solid  cylinder  with  Type  II  boundary  conditions 
and  with  a  uniform  initial  temperature.  Schneider's  solution  is  devel¬ 
oped  in  terms  of  zero  and  first  order  Bessel  Functions  of  the  first  kind. 
The  results  are  expressed  as  dimensionless  temperature  ratios  with  the 
Fourier  modulus,  a8/RQ^,  and  the  Biot  modulus,  hR0/K,  as  parameters. 

The  present  method  was  used  to  generate  comparative  results  that 
were  then  nondimens ionalized  and  compared  with  Schneider's  results. 

The  outside  surface  film  coefficient,  ho(0),  and  outside  fluid  tem¬ 
perature,  Tgo(0),  were  specified  as  arithmetic  statement  functions  in 
the  program.  The  outside  film  coefficient  was  held  constant  for  a  given 
run  but  was  varied  between  runs  so  as  to  allow  variations  in  the  dimen¬ 
sionless  Biot  modulus,  hR0/K.  A  constant  fluid  temperature  of  800°R 
was  used  in  the  calculations. 

A  comparison  of  the  present  method  with  the  exact  method  is  pre¬ 
sented  in  Table  II.  The  results  of  Schneider's  analysis  for  selected 
values  of  the  Biot  and  Fourier  moduli  are  presented  along  with  the  re¬ 
sults  of  the  present  numerical  solution.  It  is  apparent  from  Table  II 
that,  for  the  cases  studied,  a  close  degree  of  correspondence  exists 
between  the  exact  and  the  recursive  solutions. 


5.2  SAMPLE  PROBLEM  TWO 

In  Sample  Problem  Two  application  is  made  to  an  interesting  prob¬ 
lem  which  occurs  in  a  pipeline  carrying  gas  that  is  being  discharged 
from  a  pressurized  storage  vessel.  For  a  gas  stored  at  high  pressure, 
for  example,  P  =  2500  psia,  and  at  ambient  temperature  and  with  the 
pipeline  existing  initially  at  ambient  temperature,  it  is  possible  to  in¬ 
duce  severe  temperature  gradients  in  the  pipeline  because  of  the  cool¬ 
ing  effect  the  drawdown  process  has  on  the  stored  gas.  The  combined 
conditions  of  high  pressure  and  severe  thermal" gradients  could  cause 
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serious  stress  problems  within  the  pipe  wall.  Therefore,  the  unsteady 
thermal  gradients  which  would  exist  in  the  pipe  wall  become  significant. 

The  problem  of  determining  the  transient  thermal  gradients  was  in¬ 
vestigated  using  the  current  numerical  method  for  the  case  of  a  14-in.  - 
diam  steel  pipe  carrying  air  that  was  being  discharged  from  a  storage 
vessel  at  a  constant  mass  flow  rate. 

5.2.1  Program  Inputs 


The  initial  storage  pressure  and  temperature  of  the  air  was  3600 
psia  and  576°R,  respectively.  The  temperature  of  the  air  was  found  to 
be  a  function  of  the  vessel  volume,  initial  storage  temperature,  rate  of 
mass  withdrawal  from  the  tank,  and  the  time  lapse  since  initiation  of 
drawdown.  Assuming  the  air  in  the  tank  expanded  isentropically  gave  a 
relationship  for  the  air  temperature  is  a  function  of  time  as  follows: 

_  Tgi(0)  =  576[1  -  .006179  6]‘ 4  °R  (18) 


The  derivation  of  Eq.  (18)  assumes  that  the  section  of  pipeline  to  be 
investigated  is  located  immediately  downstream  of  the  storage  vessel  at 
a  point  where  the  gas  has  not  had  time  to  exchange  heat  with  the  pipeline. 
The  air  temperature  surrounding  the  pipe  was  assumed  to  remain  con¬ 
stant  at  560°R.  The  inside  film  coefficient  was  assumed  to  follow  the 
relationship  given  by  Eq.  (19). 


35 


Btu 


ft3— sec— °R 


(19) 


For  the  particular  blowdown  rate  (M  =  805  lbm/sec),  Eq.  (19)  becomes 


hi(0) 


1.6919  (rRo)-1-8  (K),6S  Qz)-'*5  (Cp)-35 


Btu 

ft3— sec— °R 


(20) 


For  the  particular  temperature  range  of  interest,  300  to  800°R,  the 
specific  heat,  viscosity,  and  thermal  conductivity  of  the  air  have  been 
found  to  obey  the  following  relationships: 

Specific  Heat 


Cn  = 


Viscosity 


/i  = 


=  0.2431  -  .00001919T  +  .0000000244T3  (  Btuc  (T  =  *11)  (21) 

V  *°m—  R/ 

(T  =  °R)  „  (22) 

W,rt* 

iv  '1 


7.3026, T3/2' 

■  x  10 


(T  +  198.6) 


sec— ft 
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Thermal  Conductivity 

i 

K  =  0.003716  +  .00001676T  +  0.000000008T’  B~ —  (T  =  °R)  (23) 

"  hr — 1 1—  R 


The  outside  film  coefficient  is  assumed  to  correspond  to  natural 
convection  from  a  horizontal  cylinder  and  is  given  by  Eq.  (24). 

ho(0)  =  0.00005  (,Tgo  -  Tn')1/3  -{t,_Bt^R"  (24) 

Equations  (18)  through  (24)  were  input  into  the  computer  program  as 
arithmetic  statement  functions.  Figure  3  shows  the  input  cards  for  the 
program  constants,  while  Fig.  4  shows  the  input  cards  for  the  initial 
temperature  distribution.  Tj  corresponds  to  the  innermost  element. 

5.2.2  Results 

The  results  of  the  study  are  shown  in  Fig.  5,  where  the  tempera¬ 
ture  as  a  function  of  location  within  the  pipe  wall  is  plotted  for  selected 
time  intervals. 


SECTION  VI 
DISCUSSION 


The  preceding  sections  have  presented  the  development  of  the  re¬ 
cursive  method,  its  application,  and  a  comparison  of  typical  results 
with  a  corresponding  exact  solution.  The  comparison  indicated  close 
agreement  between  the  recursive  method  and  the  exact  solution  for  the 
particular  cases  that  were  investigated. 

The  recursive  technique  presented  is  particularly  useful  in  cases 
involving  complicated  or  noncontinuous  boundary  conditions,  since  the 
technique  will  handle  any  boundary  condition  that  can  be  expressed  as 
either  a  function  of  time  or  a  function  of  the  surface  temperature.  The 
technique  also  allows  any  type  of  initial  temperature  distribution  to  be 
specified. 

Computer  time  required  is  not  excessive.  On  the  average,  a  nom- 
x  inal  time  of  0.  15  sec  was  required  to  compute  the  temperature  distribu- 

«•>  tion  over  a  single  time  increment  in  a  system  composed  of  ten  elements 

Since  the  recursive  technique  is  designed  to  step  through  the  solu¬ 
tion  in  discrete  time  intervals,  a  certain  inaccuracy  is  inherent  in  the 
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early  stages  of  the  solution.  These  inaccuracies  disappear  as  the  solu¬ 
tion  progresses.  In  general,  accurate  results  can  be  expected  at  any 
time,  0,  provided  the  condition  of  Eq.  (25)  is  met. 

6  l  NA0  (25) 


REFERENCES 


1.  Schneider,  P.  J.  Conduction  Heat  Transfer.  Addis  on -Wesley 

Publishing  Company,  Inc.,  Cambridge,  Massachusetts,  1955. 

2.  Eckert,  E.  R.  G.  and  Drake,  R.  M. ,  Jr.  Heat  and  Mass  Transfer. 

McGraw-Hill  Book  Company,  Inc.,  New  York,  1959. 

3.  Boelter,  L.  M.  K. ,  Cherry,  V.  H. ,  Johnson,  H.  A.,  and 

Martinelli,  R.  C.  Heat  Transfer  Notes.  McGraw-Hill  Book 
Book  Company,  Inc.,  New  York,  1965. 

4.  Smith,-  G.  D.  Numerical  Solution  of  Partial  Differential  Equations. 

Oxford  University  Press,  New  York,  1965. 


12 


AEDC-T  R-69-87 


APPENDIXES 

I.  DERIVATION  OF  THE  GOVERNING  DIFFERENTIAL  EQUATION 

II.  ILLUSTRATIONS 

III.  COMPUTER  PROGRAM  FORTRAN  LISTING 

IV.  TABLES 


13 


AEDC-TR.69-87 


APPENDIX  I 

DERIVATION  OF  THE  GOVERNING  DIFFERENTIAL  EQUATION 


Figure  6  shows  the  geometrical  aspects  of  the  mathematical  model. 
The  differential  equation  of  the  problem  is  developed  by  making  an 
energy  balance  on  a  typical  differential  volume  element  of  the  pipe  wall. 

The  following  assumptions  are  made  concerning  the  problem: 

1.  All  energy  transport  occurs  through  thermal  conduction  except 
at  the  boundaries. 

2.  Energy  transport  occurs  in  the  radial  direction  only. 

3.  Pipe  thermal  conductivity  is  constant. 

The  energy  balance  is  expressed  in  simple  terms  as 

Energy  In  —  Energy  Out  =  Energy  Stored 


The  energy  conducted  into  the  volume  element  is  given  by 
Conduction  In  =  -  ^KA  -^0  =  -  ^KRd<£dZ  -^0 

The  energy  conducted  out  of  the  volume  element  is  given  by 


Conduction  Out  =  —  (KA  — + 

\  <?R  dR 


dT  .  d 


[“  ^]"R) 

=  -  ^KRd^dZ  |KRd<^dZ  ~~~]dR^ 

=  -  KRdcidZ  M  -  Kd<£dZ  -£■  [r  -ill  dR 

<?R  <?R  L  <9rJ 

Combining  Eqs.  (1-1)  and  (1-2)  gives 

Conduction  In  -  Conduction  Out  =  KdcidZdR  IR  -  +  -^X| 

L  <?RJ  <?RJ 

The  energy  stored  is  given  by 

Energy  Stored  =  p  V  C  =  p  Rd<^dZdRC 

dd  36 


(1-1) 


(1-2) 


(1-3) 


(1-4) 
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Equating  Eqs.  (1-3)  and  (1-4)  gives  the  governing  differential  equa¬ 
tion: 

(?T  _j_  dT_ _ i_  dT_ 

dR1  +  R  3R  ~  a  36  (1“5) 


K 

where  a  =  —  and  is  known  as  the  thermal  diffusivity. 


Associated  with  Eq.  (1-5)  are  the  following  initial  and  boundary  con¬ 
ditions  : 

(1)  Initial  Condition 

at  6  =  0,  T(R,0)  =  T(R,0) 


(2)  Boundary  Condition  1  (Inside  the  Pipe) 


at  R  =  R;, 


,,  <?T(R,0) 

K-iir 


K 


<9t(r  ;.fl) 

<?R 


Here  h]{9)  and  Tgi<0)  are  the  time -dependent  inside  fluid  heat-trans 
fer  coefficient  and  inside  fluid  temperature,  respectively.  Rewriting 
this  boundary  condition  gives: 


<?T(rif  6) 
dR 


-  T(R;,0)  =  -  Tgi(0) 


( 3)  Boundary  Condition  2  (Outside  the  Pipe) 
at  R  =  R0|  K  -dT(R’  =  K.  3  T(R°’  6 


dR  dR 

=  K(9)  (Tgo(0)  -  T(R0,  6) 


Here  ho(0)  and  TgO(0)  are  the  time -dependent  outside  fluid  heat- 
transfer  coefficient  and  outside  fluid  temperature,  respectively.  Re¬ 
writing  this  boundary  condition  gives 


+  -o(^L T(R  fl) 
dR  K- 


ho(0) 

K 


Tgo(0) 


Inspection  of  the  boundary  conditions  reveals  that  they  are  not 
homogeneous;  consequently,  the  path  to  an  exact  solution  does  not 
readily  present  itself.  Recourse  must  therefore  be  made  to  approxi¬ 
mate  numerical  techniques. 
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0 ,  sec 

Fi  g.  2  Effects  of  Convergence  Criteria  on  the  Convergence  of  the  Numerical  Solution 
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Fig.  3  Sample  Problem  Two  -  Program  Constant  Inputs 


A  E  DC-T  R -69-87 


AEDC-T  R -69-87 


512.40 


0000000  oo 

i  <  >  i  ■  ■  i 

111  11  I  I  I  ill 


2222  22222 
23333  3333 
*44444  444 
55  5555555 
66GGGG66G6 
7777777777 
88010  0880 
9999995999 

i  ]  1 1  m  t  v 
GlCIt 


511.52 


0  0  0  0  0  0  0  0  0  0 

|ll  1711 II  till  11  II  If  71£I 

II  Mill 


5t  1 . 08 


000000  000 

h  n  a  a  n  a  n  a] 

III  I  I  I  I  I 


2222222  2 
3  3  3  3  3  3  3  3  3j 
444444444 
55  555  555 
GGGGGGGGGGl 
7777777777 

8  8  8  8  8  0  0  8  7] 

9  9  9  9  9  9  9  9  9  j| 

I  17  II 14 11 II  II II 


2|2  2  2  2  2  2  2  2  2  »  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  7 
3  3  3  3  3  333  313  3333  3  3  3  3  3  3  3  3  3  3  3  3  31 


■510.84 


0000  8000 


510.80 


3  0  0  0  0  0  0  0, 

jii  E  u  a  a  x  e  a  a  0U1 11 41 44 «  *  4;  a  a  a| 

111  1  I  I  1  I  It  I  1  I  I  1  1  1  I 


000000000000000000000000000008 
jsi  an  14  a  mu  si  naii  au  an  u  nan  a  11  an  h  a  an  ana 
I  1  1  I  I  1  1  1  1  1  1  1  1  1  1  1  I  1  1  1  1  I  1  1  1  1  I  1  1  1 


*|*  4  4  4  4  4  ilnitial  Temperature  i|4 
Distribution,  °R 
55  555*Format  (FIO.O) 

GGGGGGGGGGl 


7  77777777  7| 

88888  8  88 


666668  666 
7  7  77  7  7  7  7  7  ? 
00888  8888 


99999999  9|9  99999999  P|"  99999999  5 

11  aii»aaii  a  a;. |i-,  r  lima  n  a  a  caH4ian  «4i.i| 

40  1014  1011 


666666666 
7777777777 
80888  8809 


222222222222222222222222222222 
333333333333333333333333333333 
44444444444444444444444444444 
5555555555555555555555555555555 
66 66 666 6 6666 6 66 66 8686 6666 8 669 B 6 
777777777777777777777777777777 
888888880000080886860080000009 


999899999999999999999999999999 

piuuassisisianiiRBMusiiiBsiiinniinniiTiiiiia 


519. 


aoo  010080 

1  1 1  1 1  I  1  1  ■  a 

I  1  1  1  1  1  1  I  1  111 


9  9  9  9  9  9  9  9  9|9 

niiuir 

GlOM  NO 


■509. 


000  000000 
]ii  11 11 1<  n  a  11  n  a 
I  I  1  I  I  1  I  1  111 


1 1*1 


999  99999 

i  u  11  w  11  a  11 11  a  api 

11ANI 


5H3.38 


*8 

000  ooooool 

!iEn»aanana| 

11111111 


000  ooooop: 

|ii  v  »  a  a  x  v  a  a  «|n 
11111118111 


|9  99999999  919 

iiR&ana&aaMpi 
ID  iCIM  5011 


•508.  31 


■508.00 


l10 

00  0  0  0  Dll 

11  04)44 «  anuaiq 

I  1  1  1  I  I  1  1  1  1 


9999999999 

iixixiasaxa 


222222222  2|2  222222  2  2|Z  22222222  2)2  22222222  2]2  222222227222222222222222222222222222222 
33333  3  3  3  3|3  3  3  3  3  3  3  3|3  3  3  3  3  333  3)3  3333  |  333  3)3  3333  3333333333333333333333333333333333 

444444444  4|4  44444444  4|4  44444444  4I4  44444M4  4)4  444444444444444444444444444444444444444 
5  5  5  5  5  5  5  5  5)5  5  5  5  5  5  5  5  5)5  5  5  5  5  5  5  5  5(5  5  5  5  5  5  5  5  5|5  5  5  5  5  5  5  5  5)5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5 
6666666  6668366666066666666666666666666666666666661666666666666666666666668666066 
777777  7777777777777777777777777777777777777?77777|777777777777777777777777777777 
8  8  0  0  8  8  0  8  8|0  6  0  8  0  8  8  8  018  6  8  8  8  8|6  0  8  8  8  8  8  818  8  8  8  0  8  80|B  8  8  08  8  0  8  8  0  0  8  8  8  0  1  8  0  8  8  8  0  0  8  8  8  8  8  8  8 


999999999 

11  vcaaaiiaautii 


000800000000000800000D0800D600 

[iiximxxiiaxaiiDiiiiBxnxBaiinxxxnnxxx 

111111111111111111111111111111 


999999999999999999999999999999 
1  a  u  »  a  a  a  a  a  a  11  o  a  a  t>  a  n  a  n  a  11  n  » 14  n  a  »  a  a  v 


Fig.  4  Sample  Problem  Two  ■  Initial  Temperature  Distribution 
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APPENDIX  III 

COMPUTER  PROGRAM  FORTRAN  LISTING 


£***************  *******  ***  ******  **  *****  **  *  ******  *****  *****  4c 
£***  ♦  *»*»♦*«»*»*»*»»**»»♦***»♦*♦*:*  *»»**»»»»»*»  *****  ******** 
c *******  *********** *******  ********  ********* ******* ********* 

C  FORTRAN  LISTING  OF  PROGRAM  AS  APPLIED  TO  SAMPLE  PROBLEM 
C  TWO.  SEE  SECTION  IV. 

£**«»»»* »**»»*»♦«»****»»»»»♦*»**»♦ ****♦»♦♦♦ ******* *** ****** 

C* ************************************ ********************* 

C  ANALYSIS  OF  TRANSIENT  TEMPERATURE  01  STRIBUTI 0N__ IN.  A_HO L L OW 
C  CYLINDER  WITH  NON-STEADY  BOUNDARY  CONDITIONS. 

£*  ♦  ***********  ******  ******************  ************ ********* 

C  INSIDE  FILM  COEFFICIENT  ( HI >  IS  FOR  AIR-  TURBULENT  FORCED 

C  CONVECTION  WITH  VARIABLE  PROPERTIES. _ 

C **************************************************** ****** 

C  OUTSIDE  FILM  COEFFICIENT  I  HO )  IS  FOR  AIR  -NATURAL _ 

C  TURBULENT  CONVECTION  OVER  A  HORIZONTAL  CYLINDER. 

C**»**»»»***»***«»* «♦*♦♦»»♦♦*»«*♦**♦***»**»»»**♦ ♦♦***»*»»** 

C* ************************************************ ********* 

DIMENSION  A(  20  )  ,B  ( 20  ) ,  T I  TLE  (  20)  _ 

TGlJ  (  E  )  =  560.+l.*E-i.*E 

a3TGI(E)=  576. *(  (1.-.006179*E)**.4)  _jAjn 

y,jtfvVC  PGA  S  ( TG I  >  =.2431-  .000019 19*  TGI  .  0000000244*1  TG 1**2 ) 

^<AMQ'.(  TGI)=7.3026*(TGI**(3./2.  ) )  *  .  0000001/ 1  TGI  ♦198.6)  -  H 

AKGAS  ( TG I ) = (  .00  37  16+  . 00 001676* TGT+. 00 008*(TGI  **2  )  *  - „  *  * 

_ l  .«D.QP_U/3600Li _ _ _ 

- — — ^ H I  (CPGAS,AMU,AKGAS)=(  1.6919/ ((  R*RO)  **1.8  )  )  *(  AKGAS** 

1.65)*(CPGAS**.35)/t  AMU**. 45  I _ _ _ _ 

yJ>  TGOtTN)  =.00005*  (  (  ABS  (  TGO-TN)  )  **(  1  -  /  3  .  J  )  VY 

£_  F1(U.V)=(  U-V )  /  ALOGI  ( <2.*K+1.  )  +  (  2  .*N-2  .»K- 1 .  )  *R)  /  (  (2.*K 
1-1 . )  +  ( 2.*N-2.*K+1. )*R) ) 

F2  (U.V  )  =  l  U-V)  /ALOGI  (  (  2.  *K-1 .  )  +  (  2  .*N-2  .*K  +  1 .  l*R.)./.(  J  2  .  *K 
l-3.)+(2.*N-2.*K+3.)*R) ) 

F3(U)=(1.-U)*(  (2.»N-2.*K+1. ) *R  »(2.*K-1.)) 

RE  AD  (  5  «  1 )  TITLE 

1  FORMAT ( 20A4) _ 

WRITE  (6,1)  TITLE 

KEAD(5.2)BFTA.D.F,ALP.R.R0.N  _  _  _ 

2  FORMAT (6F10.0, 12) 

WRITE  (6, 3)  D.F,  ALP  .BETA, R  ,R0,  N _  _ 

3  FORMAT (/////•  INPUT  CONSTANTS'///'  DELTA  TIME  ■  'F6.3 
1'  SEC.'/'  TOTAL  RUN  TIME  =  ' F7 . 3 ' SEC . ' / '  THERMAL  DIFF 1 
2,'USIVITY  =  »F10.8'SQ.FT./SEC.'/'  BETA  =  'F10.4*  SO.', 

3 'FT .-SEC.-R/BTU'/ '  RADIUS  RATIO  (RI/RO)  =  'F5.3/  _ 

4*  OUTSIDE  RADIUS  (RO )  =  'F6.3'FT.'/'  THE  NUMBER  OF', 

5 «  ELEMENTAL  CYLINDRICAL  SHELLS  IS  (' 12' )♦///) _ 

READ (5,4) { A( J) ,J=l,N) 

4  FORMAT  (5 FI 0.0) _ 

WRITE  (6,5) 

5  FORMAT (3X« TIME '4X'T( 1 ) ' 4 X« T( 2) * 4X' T( 3 ) ' 4X* T( 4) '4X 
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1 ' T  (  b  ) '  4X ' T  (  6  )  *  4  X '  TI 7) ,4X,T(d),4X,T{9) ' 4X ' T ( 10) *  4X 
.I'LGAS.I  '.WTGASO’  ) 

CHI=2.*(N**2)*D*ALP/  ( RJ*  *2 ) 

_ E  =  0 . 0  _ _ 

b  T  G  A  S  0  =  T  G  0  (  F  }  .  . . . 

TGAS I=TGl{ E  ) 

WRITE  (6,7)  E,C  A(  j),J  =  L,N),  TGASI  ,  IGA  SO 
7  FOE  M  AT  ( /  l  3F_8  . 2  ) 

'  R  a  i" 

_ b  (JJ  =  CHI  *  (  (R*BEjn\*H[(CIJGAS{  FGI  (  E  ) )  ,AMU(TGI  (6  ) )  » AKGAS  ( 

1TG  I  (E)  )  )  *  (  TGI  ( E  )-a'(T>  )  T+F  1C  M2),  A(  f)  )  )  /F3t  R  )  +  All) 
NN»N-1 
00  B  K=2  t  NN 

d(K)=  CHl*(F2( A(K-L) ,A{K))+  FI ( A (K+ 1 ) , A( K> )} /F3 ( R ) 
UA(K) 

3  CONTINUE . . .  . . 

K  =  N 

B(N)=CHI*( F2IA IN-1) ,A(N) ) +B2  TA*HO( TGO ( E ) ,A{N) )*(TGD(E) 
L-AIN)  ))/F3(R)+A(N) 

00  9  J=1 ,  N 
A  (  J ) =B ( J ) 

9  CONT  I_NUE _ _ 

E  *  E  +  0  . 

I F ( E-F 16,6,10 
10  STOP 
ENH- 
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TABLE  I 

ERROR  PROPAGATION  IN  A  TWENTY-ELEMENT  NUMERICAL  SYSTEM 


Element  No. 


Time  (sec) 
0.  0 
0.001 
0.  002 
04 
08 
15 
30 
45 
0.  060 
0.080 
0.  100 


0.02 

0.03 

0.04 

0.07 

0.  07 

0.07 

0.07 

0.07 

0.07 

0.06 

0.06 

0.06 

0.04 

0.  04 

0.04 

0.03 

0.03 

0.  03 

0.  02 

0.  02 

0.  02 

0.02 

0.02 

0.02 

0.01 

0.01 

0.  01 

001  sec 

0  0.01  0.02 

0. 04  0.  05  0.  07 

0.07  0.08  0.08 

0.07  0.07  0.07 

0.06  0.06  0.06 
0.04  0.04  0.04 

0.03  0.03  0.03 

0.02  0.  02  0.02 
0.02  0.02  0.02 
0.01  0.01  0.01 


11 

12 

13 

0 

0 

0 

0.  15 

0. 12 

0.06 

0.  10 

0.  0!) 

0.  07 

0.  07 

0.07 

0.06 

0.  06 

o.  or. 

0.05 

0.  05 

0.  04 

0.04 

0.  04 

0.03 

0.  03 

0.  03 

0.  03 

0.  03 

0.  02 

0.02 

0.  02 

0.  02 

0.02 

0.03 

0.  01 

0.01 

0.01 

0 

0.  01 


16 

17 

18 

19 

20 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0.  02 

0.  01 

0.01 

0 

0 

0.  03 

0.  03 

0.02 

0.02 

0.  02 

0.04 

0.  03 

0.03 

0.  03 

0.03 

0.04 

0.  04 

0.03 

0.  03 

0.03 

0.  03 

0.  03 

0.  03 

0.03 

0.03 

0.  02 

0.  02 

0.02 

0.02 

0.02 

0.  02 

0.  02 

0.02 

0.  02 

0.02 

0.01 

0.  01 

0.01 

0.01 

0.01 

0.  01 

0.  01 

0.01 

0.  01 

0.01 
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TABLE  II 

COMPARISON  OF  RECURSIVE  SOLUTION  WITH  EXACT  SOLUTION 
FOR  TRANSIENT  TEMPERATURE  DISTRIBUTION  IN  A  SOLID  CYLINDER 


Biot  Modulus, 
hR0 

K 

Fourier  Modulus, 
ad 

R02 

Radius  Ratio, 

R  /  R(j 

Dimensionless 

Temperature  Ratio, 

T/ ^initial 

Exact  Solution 

Recursive  Solution 

0.  1 

0.  1 

0 

1.  00088 

1.  00085 

0.  1 

0.  1 

1 

1. 01352 

1.  012767 

0.  1 

1.0 

0 

1.  05234 

1.05227 

0.  1 

1.0 

1 

1.  06588 

1.06518 

1.0 

0.  1 

0 

1. 00798 

1.  00803 

1.0 

0.  1 

1 

1. 10515 

1.  10202 

1.0 

1.0 

0 

1.  25021 

1.  25*238 

1.0 

1.  0 

1 

1.  27989 

1.  28050 

10.  0 

0.  1 

0 

1.  03331 

1.  03618 

10.  0 

0.  1 

1 

1. 28946 

1.  28987 

10.0 

1.  0 

0 

1.  32881 

1.  32967 

10.0 

1.0 

1 

1.  33278 

1. 33288 
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